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ABSTRACT 

The analysis of the vertical velocity dispersion of disc stars in the local Milky Way is the most direct astronomical means of estimating 
the local dark matter density, poM- Current estimates for poM based on the mid-plane dynamic density use a local baryonic correction 
that ignores the non-local effects of spiral structure and significantly underestimates the amount of dynamically relevant gas now 
known to be present in the ISM; the additional gas plus the remaining uncertainties make it practically impossible to measure poM 
from mid-plane kinematics alone. The sampling of inhomogeneous tracer populations with different scale-heights and scale-lengths 
results in a systematic increase in the observed dispersion gradients and changes in the nominal density distributions that, if not 
properly considered, can be misinterpreted as a sign of more dark matter. If the disc gravity is modelled locally using an infinite disc, 
the local variation in the vertical gravity due to the globally exponential disc components results in an underestimation of the baryonic 
contribution by as much as ~40%. Given only the assumptions of stationarity, an axially and vertically symmetric disc, doubly 
exponential tracer and mass-component density profiles, a phenomenologically justified model for the cross-dispersion component 
ctr;, and a realistic model for g-iz), it is possible to solve the full vertical Jeans equation analytically for the vertical dispersion cr^fz) 
and hence test the robustness of previous attempts at measuring poM- When the model parameters for are estimated from SEGUE 
G dwarf star data, it is still not possible to explain the difference in behaviour seen in the simple thick- and thin-disc datasets reported 
by Buedenbender et al. (2014). Rather than being a fundamental problem with the kinematical model, this effect appears to be a 
further sign of the difficulty of defining and handling kinematically homogeneous tracer populations. 
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1. Introduction 


Although there are many indirect lines of evidence for the pres¬ 
ence of non-baryonic matter in the Universe, nowhere is a con¬ 
crete estimate for the density of such dark matter (DM) more im¬ 
portant than in the Solar neighbourhood: the local density deter¬ 
mines how easily such matter might be detected in laboratory ex¬ 
periments. Starting with the classic analyses of Kapteyn (1922) 
and Oort (1932, 1960), there have been many attempts to use the 
kinematics of stars in the local Milky Way (MW) to estimate 
the local density and distribution of mass. Recently, the data 
required for such an analysis have become much better thanks 
to massive kinematic surveys like SEGUE (Yanny et al. 2009), 
RAVE (Steinmetz et al. 2006), and APOGEE (Allende Prieto 
et al. 2008). Read (2014) has recently given a detailed review 
of the theoretical expectations from a cosmological/extragalactic 
perspective and the history of the astronomical attempts to mea¬ 
sure the local DM mass density; he shows that recent estimates 
have yielded values in the range puM = 5 - 15mMopc“^, but 
that systematic effects are still the dominant source of error|^ 
The measurement can be made using one of two different basic 
approaches, each with its own strengths and weaknesses. 

One can attempt to fit the global kinematics of the stars and 
gas in the MW in detail and, by generating a global mass model, 
also estimate the local DM density (e.g. Bienayme, Robin & 
Creze 1987, McMillan 2011). The problem with this approach 
is that basic properties like the distances to kinematic structures 
(e.g. Binney & Merrifield 1998 Section 9.1.1; Reid et al. 2014 ) 


* I use the astronomical unit 1 mMopc^^ = 0.001 Mopc^^ = 

0.038 GeV to avoid having to count the many zeros present when 
using Mopc“^. 


or even the distance to the center of the MW (e.g. Gillessen et 
al. 2009; Sofue et al. 2011; Branham 2014 ) are difficult to mea¬ 
sure; the area within which a reasonably robust analysis can be 
made is only between the end of the central bar and the Sun’s 
Galactocentric radius, Ro (Dehnen & Binney 1998); and the pa¬ 
rameters of arbitrary models for the DM geometry must be fit, 
e.g. using a Navarrow-Erenek-White (1997) or Einasto (1965 ) 
profile plus assumptions about the halo’s oblateness. A recent 
example of this approach is that of Piffl et al. (2014): using 
~ 200,000 stars from the RAVE survey with Galactic heights 
|z| < 1.5 kpc and a 33-parameter model for the global mass- 
distribution and its effects on the stellar kinematics via assumed 
stellar distribution functions (DE), they derived a local DM den¬ 
sity of 15 + 2mMopc“^. When kinematic data from the GAIA 
satellite become available, this method will be the only way of 
dealing with the complexity of a global data sample from a re¬ 
alistic MW (Binney 1998 ) and the sheer enormity of the dataset 
will constrain the otherwise arbitrary model parameters to con¬ 
form to reality. 

An alternate approach better suited to the limits of current 
data and our knowledge of the structure of the MW (Carraro 
2015 ) is to try to measure the DM density locally using local 
stellar kinematics (e.g. Bahcall 1984 Kuijken & Gilmore 1989c 
Holmberg & Elynn 2004 Bovy & Rix 2013 ). The advantage of 
this approach is that the proper motions and distances to all stars 
and hence the full three-dimensional densities, velocities, and 
dispersions can principally be measured, at least at large dis¬ 
tances from the plane of the MW where the effects of DM should 
show up most directly, and the large-scale distribution of DM 
presumedly plays a minor role (e.g. Siebeit et al. 2008). Here too 
the complex potential and DE models can be used, but they may 
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not be able to ensure that the result is accurate (versus merely 
precise). For example, the simple two-component Stackel poten¬ 
tials used by Bienyame et al. (2014) to model the disc and DM 
halo give the local disc a power-law rather than an exponential 
radial profile; any inaccuracies in the resulting model dispersion 
will be compensated by modifying the ill-constrained properties 
of the DM halo. 

It is important to try to estimate the gravitational effects at¬ 
tributed to Pom locally on length-scales smaller than Rq, since 
there are competing theories like MOND (e.g. Milgrom 2014) 
that invoke changes in the behaviour of gravity on scales of 
whole galaxies - the same scales relevant for studies of the 
MW’s rotation curve. If the local effects of DM in the MW can 
only be determined using the MW’s rotation curve, we will have 
lost the information at one scale-length and have to deal with 
much more complexity and intrinsic uncertainties. The purpose 
of this paper is then to probe how robust truly local estimates 
for Pom can be. I first show that the extent to which the derived 
values of poM are compromised by the many uncertainties and 
pitfalls in the methods and data used is generally not sufficently 
recognized. Then, I describe a relatively simple method with 
which one can principally make more robust estimates for the 
local DM density using the full vertical Jeans equation; as few 
assumptions as possible are made about the large-scale stellar 
kinematics, MW properties, and unknown DM model parame¬ 
ters, and the effects of radial coupling and a realistic rather than 
an idealized toy model is used for the local gravity. Finally, I ap¬ 
ply this method to a particular dataset in order to see whether the 
challenges fundamentally lie in our ability to model the physics 
or to have adequate data. As we shall see, the problems with 
many previous poM estimates are serious and one must more 
carefully consider the inadequacies of models and the difficulty 
of defining and handling appropriate tracer datasets. 

2. Near-midplane estimates for the DM density 

If the stellar motions and the local potential <l)(z) near the 
MW mid-plane are simple enough, the Jeans theorem says that 
Vj/ 2 H-<t) is an approximate isolating integral and the dynamic lo¬ 
cal density p^y„( 0 ) can be calculated from the distribution of mid¬ 
plane vertical velocities V 2 ( 0 ) and the observed density distribu¬ 
tion v(z) of a tracer population (Kuijken & Gilmore 1989c; see 
Sect. 4.9.3 in Binney & Tremaine 2008 hereafter BT08). Given 
enough good data, this method should yield a very robust esti¬ 
mate for pdyniO) and - given a sufficiently accurate model for the 
baryonic contributions - a similarly robust estimate for the re¬ 
maining Pom- This analysis has been made using many datasets 
and models (e.g. Kuijken & Gilmore 1989b; Holmberg & Flynn 
2000; Creze et al. 1998; Garbari et al. 2012). 

The problem with the conversion of the dynamic density, 
once measured, to the local DM density is the difficulty of es¬ 
timating the baryonic contribution pdyn,bi^)- The observed mid¬ 
plane density of stars is fairly simple to estimate from local 
surveys; Read’s (2014) re-tabulation of the Flynn et al. (2006) 
model has 37.9 and 3.5 mMopc^^ from the thin and thick discs, 
respectively (if one assigns the local giants, white dwarfs, and 
brown dwarfs to the thin disc for this purpose). If one assumes 
that the density-ratio of the thick disc to thin disc is the value 
measured photometrically by Juric et al. (2008), then the contri¬ 
butions have to be redistributed as 37.0 and 4.4mMopc“^. The 
local surface density of molecular, neutral, and warm gas has 
been classically estimated as 13Mopc“^ by Holmberg & Flynn 
( 2000 ), from which one can estimate the local densities using 
typical observed vertical scale-heights and/or linewidths. The 


Table 1. A local baryonic mass model using classic (updated) ISM den¬ 
sity estimates. 


Description 

P(0)'' 

O'! 

^HF* 


(mMopc^^) 

(kms-') 

(Mopc-2) 

Molecular gas 

21.0(35.0) 

4.0 

3.1 (5.1) 

Cold neutral medium 

16.0 (32.0) 

6.1 

3.6 (7.3) 

Warm neutral medium 

12.0 

14.8 

7.3 

Warmer gas 

0.9 

40.0 

2.0 

Mv < 2.5 

3.0 

7.5 

0.8 

2.5 <Mv <3 

1.5 

10.5 

0.6 

3 < Mv < 4 

1.9 

14.0 

1.1 

4 < Mv < 5 

2.1 

18.0 

1.6 

5 < Mv < 8 

6.8 

18.5 

5.4 

My > 8 

13.1 

18.5 

10.5 

Giants 

0.6 

20.0 

0.5 

White dwarfs 

5.8 

20.0 

5.1 

Brown dwarfs 

1.9 

20.0 

1.7 

Thick disc 

4.4 

37.0 

8.6 

Stellar halo 

0.1 

100.0 

0.7 

Sum of gas 

50(80) + 25 

- 

16 ( 22 ) 

Sum of stars 

41+4 

- 

37 

Total sum 

91(121)+ 25 

- 

53 (58) 

’ Based on Flynn et al. (2006); see text for details. 


*Calculated assuming isothermal components and the potential 

of Holmberg & Flynn (2004 ). 



stellar densities are generally thought to be good to about 10 %, 

whereas the gas densities 

are difficult to estimate; the latter are 

given errors of ~50% by Homberg & Flynn (2000 ), suggesting a 
total error of about 25 mMopc^^. It should be noted that the lo- 


cal baryonic mass-density within this traditional model is dom¬ 
inated by the gas - i.e. by the least well-determined component 
- and the uncertainty is already much larger than the remaining 
5-15 mMopc^^ thought to be due to DM. 

The total estimated baryonic density pb{0) ^ 91mMopc“^ 
can be converted to total surface density by adding up the contri¬ 
butions of individual kinematically isothermal components and 
hence depends upon the assumed or fitted potential. Using the 
potential from Holmberg & Flynn (2004), one obtains the val¬ 
ues shown in Table 1. Like Read, I have used updated HI disper¬ 
sions (Kalberla & Dedes 2008). These surface density estimates 
are not internally-consistent - one would have to repeat the fit 
of the K giant stars performed by Holmberg & Flynn using the 
modified densities and dispersions - but they are similar to those 
tabulated by Flynn et al. or Read and are good enough for our 
present purposes. 

The traditional estimates of the neutral and molecular 
gaseous densities from the 1980s and 1990s are off for a variety 
of reasons; we now know that the density and ionization struc¬ 
ture of the real interstellar medium (ISM) is much more com¬ 
plex. For example, a significant fraction of the HI cold neutral 
medium (CNM) is not optically thin, increasing the true densi¬ 
ties by at least ~30% (Grenier, Casandjian & Terrier 2005) and 
probably by as much as factors of 2-3 (Braun 2012; Fukui et al. 
2014). There is “dark” molecular gas not easily seen in CO but 
visible in gamma rays (e.g. Ackermann et al. 2012) and [CII] 
(e.g. Pineda et al. 2013 ) representing as much as 40% of the 
total, i.e. an increase from the traditional value of 67%. These 
effects must raise the estimates for the local ISM densities con¬ 
siderably: assuming a factor of 2 correction for the CNM and a 
factor of 1.67 for the molecular gas, Pgas,hcai(0) increases from 
50 to 8 OniM 0 pc“^. The Sun also sits near the middle of the Lo- 
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cal Bubble (e.g. Lallemont et al. 2003), a significant hole in the 
local Galactic ISM with a size comparable to the local thin-disc 
scale-height; the local gas density on somewhat larger scales is 
much higher by definition. 

Unfortunately, the derived local value of p;,(0) is not even 
what we really need: pb,dyn(0) and pbjocaiiO) are not the same, 
given that the potential of the MW disc experienced by the stars 
is highly non-uniform and that the vertical epicylic frequencies 
are ~2-3 times larger than the orbital frequency (e.g. BT08, p. 
167), effectively averaging out the local potential differences. 
The Sun sits in a trough between two major spiral arms (Perseus 
and Scutum-Centaurus) and at the edge of the Orion Spur (Xu 
et al. 2013): the local volume- and surface-densities must then 
be uncharacteristically lower than the mean kinematical values - 
again by definition - unless the effects of the Orion Spur are un¬ 
characteristically large. It is extremely difficult to measure em¬ 
pirically or estimate theoretically how big this effect is, so it is 
generally not taken into account in the literature, not even when 
fitting local or global DM models. Drimmel & Spergel (2001) 
tried to measure the spiral amplitude of the MW by modelling 
projected COBE NIR data: they placed the Sun in an inter-arm 
region with a surface density contrast '^mterarm l^arm ~ 0.85 in 
the K-band, corresponding to a spiral amplitude A - 0.14. Ex¬ 
ternal galaxies may give us a more reliable estimate for the typ¬ 
ical arm-interarm contrast of large spiral galaxies like the MW 
(Rix & Zaritsky 1994). Grosbpl, Patsis & Pompei (2004) anal¬ 
ysed ground-based K-band images of 54 normal spiral galaxies 
and found amplitudes between 0.12 and 0.23. Elmegreen et al. 
(2011) and Kendall, Clarke & Kennicutt (2014) studied 62 MW- 
like galaxies at 3.6pm using the Spitzer satellite (Sheth et al. 
2010; Kennicutt et al. 2003) and found typically larger ratios of 
A = 0.1 - 0.5. Thus, a galaxy like the MW is likely to have 
A ss 0.2, representing an interarm-to-mean surface density cor¬ 
rection of 1/(1 -A)=s 1.25. Vertical compression make the effects 
in volume density even larger than the ones in surface density: 
Dremmel & Spergel (2001) invoke a mid-plane arm-interarm 
fiux density contrast that is 1.3 times larger than the projected 
value, yielding an interarm-to-mean volume density correction 
of 1/(1 - 1.3A)=sl.35. 

The effects of spiral arms are even stronger in the ISM: while 
Holwerda et al. (2005 ) find that there is only a difference of 
2 in the opacities of arm and inter-arm regions, corresponding 
to the A ^ 0.2 seen in external stellar discs, they attribute this 
effect to cloud dumpiness and not to total surface density. In¬ 
deed, the HI spiral amplitudes seen in the THINGS datasets 
(Walter et al. 2008) are much larger than those seen for stars, 
reaching up to nearly 100% amplitudes. Thus, it is not surpris¬ 
ing that the mean HI surface densities at the Galactic radius of 
the Sun estimated by Kalberla & Dedes (2008 ) are larger than 
those of the local estimate^ Y.{HI+He) ^ 15Mopc“^ (uncor¬ 
rected for optical depth effects), consistent with a spiral ampli¬ 
tude A (15-11)/15 = 0.3._ 

Assuming thatpdy„(0) a; p(0) p/oca;(0)/(l - 1.3A), a spiral 

amplitude of 20% alone increases the effective value of pb,dyn{0) 
from 91 to 123 rriMopc^^ and with more modern estimates of the 
densities to 163mMopc“^. In order to compare these values of 
PdyniO) with those derived from recent kinematical studies (e.g. 
the exhaustive compilation in Read 2014 ), one can calculate the 
effective value of poM.eff = pdyniO) - Pb.miniO) assuming Read’s 
baryonic minimum estimate of 91 mMopc"^: either the ISM and 
spiral amplitude corrections applied alone correspond to an ap- 

^ Read (2014) forgot to correct the Kalberla & Dedes value of ~ 
12Mopc“'^ from the pure hydrogen to the total density. 


parent DM density of ~ 30mMopc“^, and both corrections ap¬ 
plied simultaneously result in ~ 72niM0pc“^. Assuming that 
the estimates for pdyniO) are accurate, it would appear that any 
reasonable correction to pi(0) for known effects is much more 
than is needed and corresponds to apparent DM densities that 
are much larger than the range 5-15 mMopc^^ usually derived 
and thought to be consistent with cosmological expectations and 
galaxy-evolution scenarios. This is already true for the tradi¬ 
tional estimate of the local baryonic density, given the uncer¬ 
tainties in the traditional estimate of the gas contribution, and the 
realistic corrections to the total ISM densities and/or corrections 
for spiral structure only make the situation much worse. Thus, 
practically any empirical value of pdyni^) in the range found by 
current kinematic studies is incapable of constraining the local 
DM density without additional (presently unavailable) informa¬ 
tion, assumptions, and/or additional external constraints. 


3. Inhomogenous data and the effects of model 
assumptions 

Since local estimates of p^M using current data are unable to 
yield robust values by themselves, we must turn to kinemati¬ 
cal effects farther away from the plane of the MW (e.g. Bahcall 
1984; Garbari et al. 2011). By far the simplest method is to use 
the Jeans equations (moments of the Boltzmann equation; e.g. 
Sect. 4.8 in BT08), which connect the observed dispersions - 
particularly the z-z component cr^(z) = cr^(z)^ - with the lo¬ 
cal gravity above and below the plane of the Galaxy. The latter 
is dominated by the relatively well-constrained stellar popula¬ 
tions, small amounts of gas, and (presumedly) a roughly constant 
DM density. Here, there are also various approaches using dif¬ 
ferent assumptions and amounts of data (e.g. Kuijken & Gilmore 
1989a; Bovy & Tremaine 2012 hereafter BT12). 

Moni Bidin et al. (2012b ) analysed the kinematic data for a 
sample of colour-selected giants taken from Moni Bidin, Car- 
raro & Mendez (2012a; hereafter MBCM), for which the com¬ 
plete dispersion tensor could be estimated, using the Jeans and 
Poisson equations to derive the surface density with height, Z(z). 
Assuming a double exponential tracer density distribution and a 
mean azimuthal velocity that is not a function of Galactocentric 
distance at any height from the plane, they derived a local DM 
density from the asymptotic behaviour of Z(z) of 0 ± 1 mMopc^^. 
Bovy & Tremaine (2012) argued that one should assume a par¬ 
ticular form for the vertical variation in the radial force field 
and a different radial scale-length, and so derived a value of 
8 + 3 mMgpc^^ from a simple vertical Jeans equation (but see 
Moni Bidin et al. 2015). Beyond difficult technical questions 
about the measurability of dynamical quantities and their gra¬ 
dients and - in the case of the full Jeans equation - the values of 
various global MW parameters, all of these analyses suffer from 
a major defect: they incorrectly assume that the MBCM data 
represent the kinematics and density of a uniform tracer popula¬ 
tion with a well-defined uniform scale-length and scale-height. 
In fact, the giants are a vertically differentiated mixture of popu¬ 
lations, each with its own very different abundance, scale-length, 
scale-height, and dispersion, and selection by colour alone is not 
enough to insure homogeneity of the spatially resolved kinemat¬ 
ics. Since the vertical dispersion is a strong function of scale- 
height, the varying mixture of tracer populations automatically 
creates a large positive gradient in the mean cr^ as well as a de¬ 
crease in the effective scale-height over the thick-disc value. Be¬ 
cause the Jeans equations are linear, spatial and kinematic in¬ 
homogeneity is not per se a problem as long as the additional 
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gradient in kinematics created by the mixture of tracer popula¬ 
tions is accompanied by an equivalent change in the combined 
density. Thus, the impact of kinematic inhomogeneity depends 
upon the model used. 

Equation 25 in BT12 - derived assuming strictly exponen¬ 
tial density distributions with scale-heights h - can be used to 
understand the effects in the MBCM data at heights |z| >> h: as¬ 
suming the change in the derived value of E(z) is mainly due to 
DM, that the effects of cross-dispersion terms are small, and that 
- consistent with the current quality of the kinematic data - any 
gradient in cr^(z) is roughly constant, one can see that the DM 
reveals itself mostly via an induced vertical gradient in cr^(z) at 
large heights; 



1 ~ M 


[ dz j 


1 9z ) 


However, if this equation is applied to an inhomogeneous mix¬ 
ture of tracer populations rather than a homogeneous thick-disc 
tracer, an increase in the gradient of a-^{z) appears due to the ad¬ 
dition of tracers with lower cr^ and smaller h. If the density of this 
mixture is fit with a single exponential vertical density profile, 
the result will be a shorter apparent scale-height and the fitted 
DM density will be higher because of the increase in the appar¬ 
ent dispersion gradient and because of the decrease in apparent 
scale-height. One can easily estimate the magnitude of this effect 
using the simple thin-Hthick disc photometric model of Juric et 
al. (2008); for two exponential components with scale-heights of 
300 and 900 pc, a local thick-to-thin density ratio of 0.12, reason¬ 
able isothermal dispersions of 15 and 45 km s ' (Kordopatis et 
al. 2013), and an observed value simply weighted by the relative 
densities, the tracer mixture effect alone is expected to produce 
a 20% increase in cr^ over the range of 1.4 to 4.5 kpc covered 
by the MBCM data, exactly as is observed. The corresponding 
change in effective scale-height depends upon the details of the 
sampling and fitting, but simple mock datasets created using the 
above model show a significant decrease in the effective vertical 
scale-height from the thick-disc value assumed by MB 12. Bovy 
& Tremaine (2012) attempted to correct for the uncertainty in 
the mean scale-height, but if an inhomogeneous distribution is 
non-exponential, the fundamental equation and not the parame¬ 
ter is the problem. Sanders (2012) has performed a much more 
detailed simulation of the MBCM dataset and analysis using DF, 
a simple thin and thick disc, an assumed evolutionary history, a 
larger (20%) local thick disc fraction, and a standard DM halo. 
He - of course - found the same increase in the gradient in cr^(z), 
but the complexity of his model has led to the impression that the 
increased gradient implies more rather than less DM (Garbari et 
al. 2012;Piffl et al. 2014). 

A similar effect plagues the analysis of Bienayme et al. 
(2014), who used RAVE observations of red clump giants, a sim¬ 
ple potential model, and an assumed DF to derive the vertical 
gravity g^{z). Fitting g^{z) a posterori with a single doubly expo¬ 
nential disc, an intermediate radial scale-length of 2.5 kpc, and a 
constant DM density, they then found poM = 14+1 mMopc^^. 
Bienayme et al. did split their kinematic dataset into three abun¬ 
dance groups, corresponding to [Fe/H] < -0.58, [Fe/H] > -0.25, 
and an intermediate population; while not enough to isolate ho¬ 
mogenous tracers, this should result in similar vertical scale- 
heights (see Fig. 4 in Bovy, Rix, Fiu et al. 2012b), so that the 
density-weighted ctj. should not be a strong function of height 
even if each subsample is kinematically inhomogenous. How¬ 
ever, interstellar extinction makes the radial range over which 
these samples were obtained a strong function of galactic lati¬ 


tude: the thin-disc data was obtained over a significantly smaller 
range of Galactic radii than the thick-disc data. Chen et al. (2012) 
showed that the thin- and thick-disc scale-lengths estimated from 
SEGUE stars near the Galactic plane are different, confirming 
the results of Bensby et al. (2011) using many fewer stars: their 
values are H,hin -3.4 and H,hick - 1.8 kpc, suggesting that the 
thin-disc scale-length is nearly twice that of the thick disc. As¬ 
suming that the squared velocity dispersions vary exponentially 
with the same scale-length and hence are larger at smaller Galac¬ 
tic radii, the thick-disc data is weighted to smaller average radii 
and hence larger dispersions. Thus, there must be an inverse cor¬ 
relation between the radial scale-length and cr^, producing an 
artificially larger dispersion and hence larger p^M for the data 
thought to be the most sensitive. Again, a crude model can be 
used to estimate the effects: for an object selection criterion 
roughly corresponding to that used by Bienayme et al., consist¬ 
ing of a minimum distance of 200 pc, a maximum distance of 
3000pc, \b\ > 22°, and radial scale-lengths of 1.6 or 3.4kpc for 
the thick and thin disc, respectively, one expects an artificial in¬ 
crease in the relative vertically averaged values of the thick-disc 
cTj over those of the thin disc over the range of 200 to 2000 pc of 
~ 20% independent of the common vertical scale-height. This 
is roughly the gradient seen in the RAVE data and a sampling 
effect not considered in their model. 

Even when considerable effort is put into an attempt to cor¬ 
rect for systematic effects, to treat consistently inhomogenous 
tracers, and to minimized the number of assumptions, the result 
is not necessarily more reliable. For instance, Garbari, Read & 
Fake (2011) use a detailed global MW model in an attempt to 
correct for the effects of a realistic Galaxy (e.g. spiral arms), 
but their mock-data simulation has neither a thick-disc com¬ 
ponent nor an ISM component even though these two compo¬ 
nents make up 60-70% of the mid-plane density (Table 1) and 
differences in radial scale-lengths between the thin and thick 
discs must also affect the sampling. In addition, their fits were 
constrained to have pbjocaii^) in the unnecessarily narrow range 
91 ± ITmMopc^^, which implies a considerable shift in kine¬ 
matic influence to the otherwise less well-constrained DM halo. 

The problems with mixed tracer populations can be mini¬ 
mized by carefully choosing each sample to be as homogenous 
as possible, and the most obvious criterion of homogeneity is de¬ 
tailed abundance, e.g. [Fe/H] and/or [a/Fe]. Zhang et al. (2013; 
hereafter Z13) derived a DM density of 6.5 + 2.3 mMopc^^ us¬ 
ing three different [a/Fe]-selected, K-dwarf tracer populations 
taken from the SEGUE survey at vertical heights between 300 
and 1400 pc, roughly corresponding to thin-, intermediate-, and 
thick-disc tracer populations. However, there are several prob¬ 
lems with this analysis as well; Z13 used an approximate vertical 
Jeans equation; they assumed that the baryonic contributions to 
the local vertical gravitational held are due to infinite homoge¬ 
nous discs; their stellar disc had a thin-disc scale-height even 
though at the higher heights the gravity should be dominated by 
the thick disc; they assumed the traditionally low gas surface 
density (see discussion in the previous section); the three (as¬ 
sumed exponential) tracer populations each cover a wide range 
in [Fe/H] and hence in scale-heights, in all probability produc¬ 
ing the kinematic inconsistencies implied by Eq.[^(the thin and 
intermediate tracers have density profiles with non-exponential 
tails, as expected); and they used an undocumented Markov 
chain Monte Carlo (MCMC) prior on the scale-heights (Zhang, 
private communication) without which their value of poM is un¬ 
determined. 

Bovy & Rix (2013 ) analysed a SEGUE/SDSS dataset of G 
dwarfs at Galactic radii 4.5 <R<1 kpc and much larger heights 
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than was possible with the K dwarfs used by Z13. They modelled 
the data using DF models for each mono-abundance ([Fe/H] and 
[a/Fe]) data subset and modified Stackel potentials, permitting 
them to simultaneously estimate important quantities like the 
radial scale-lengths. Many publications now use their derived 
value of the stellar surface density as an assumed accurate mea¬ 
sure of the total baryonic contribution (e.g. locco et al. 2015). 
However, the use of Stackel potentials results in vertical disper¬ 
sion profiles cr^{z) that increase towards the disc mid-plane for 
heights below 2kpc, and are otherwise flat (their Fig.4d); mea¬ 
surements of the near-mid-plane dispersions are significantly 
lower (Kuijken & Gilmore 198%; Creze et al. 1998; Dehnen & 
Binney 1998; Holmberg, Nordstrom & Andersen 2007; Pasetto 
et al. 2012a) and the mono-abundance data of Biidenbender et 
al. (2014; hereafter B14) clearly show that the gradient at such 
heights, if present, has the opposite sign. Part of this problem 
may be due to the fact that their model did not account for the 
effects of the local inter-arm region (see the discussion in the 
previous section), which would produce an apparent change in 
the radial scale-length and hence potentially observable changes 
in the modelled kinematics. They also used the same tradition¬ 
ally low estimate for the gaseous surface density and a single¬ 
component disc potential. Other issues are discussed by Read 
(2014). Thus, while Bovy & Rix (2013) have shown how poten¬ 
tially powerful the complicated method for modelling the kine¬ 
matics of disc stars can be, their impressively precise results are, 
in fact, inaccurate at a level that is important when trying to mea¬ 
sure Pdm- 

Biidenbender et al. (2014) separated the SEGUE G dwarf 
data into clearly thin- and thick-disc tracers and showed that the 
simple model used by K13 is unable to simultaneously explain 
the differences in the magnitudes and the slopes of the two trac¬ 
ers. They attributed this to different orientations of the velocity 
ellipsoid with height, implying that the vertical motions are not 
fully decoupled from the radial ones. However, both the thin- 
and thick-disc tracer populations were not mono-abundance se¬ 
lected and, like the Z13 samples, contained stars with very differ¬ 
ent [Ee/H] and hence vertical scale-heights. Was the difference 
due to model assumptions and fundamental kinematic effects, or 
simply due to inhomogeneous tracers? 

4. The not-so-local gravity field 

In order to derive p^M from homogenous velocity dispersion 
data using a minimally complicated model of the MW - station¬ 
ary, axisymmetric, vertically symmetric - in a manner similar to 
that used by Z13, one must specify the gravitational forces on 
the stars more accurately than possible in the simple model they 
used. Although the difference between the gravity from an infi¬ 
nite, uniform disc and the real MW should be small for |z| «Rq, 
the MBCM data reach to ~ Roll and a correction to the model 
of an infinite uniform disc that takes structure at larger scales 
is needed. Unfortunately, the effort to construct a global mass 
model is exactly what one wanted to avoid by solving the vertical 
Jeans equation alone, so the question is whether the non-uniform 
effects are significant and, if so, how easily they can be included 
as corrections to a local analysis. 

The effects due to the global variation of the disc surface 
density with radius, parameterised by the radial scale-length H 
of some component, are straightforward to calculate and are of¬ 
ten - but by no means always - included in kinematical analyses 
when they should be. Eor large H, the disc regions at smaller 
and larger Galactic radii will have similar surface densities and 
one would naively expect the correction to be small. Eor small 
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z fl^pcl 

Fig. 1. The vertical gravity of axisymmetric doubly exponential 
MW discs and radial scale-lengths H = 1.4,1.8,2.2,2.6,3.0,3.4, and 
3.8 kpc; the blue fine is for a uniform infinite plane with the default lo¬ 
cal surface density Eq = 2i(Ro; 2 I = 0); the red lines are fits to g.(z < 3 kpc) 
assuming an additional effective infinite uniform disc component (see 
text). 

H, the surface density at smaller Galactic radii will have a sig¬ 
nificantly larger surface density that could dominate the local 
vertical gravity contribution and which may or may not be com¬ 
pensated by the loss of surface density at larger radii - the grav¬ 
itational field of a planar annulus is very different from that of a 
point source, making it hard to know the outcome a priori. Eor 
a doubly exponential disc component, this effect can be calcu¬ 
lated semi-analytically using the expressions for g^(R,z) derived 
by Kuijken & Gilmore (1989c 

The results for a doubly exponential disc with vertical scale- 
height h = 300 pc and a range of H are shown in Eig. 1 (top); 

^ The vertical acceleration g, is often labelled Kp, the latter often has an 
ill-defined sign, and is often incorrectly called a “force” in the literature. 
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for H^3 kpc, is very similar to that due to an infinite uniform 
disc; for larger H there is a slight decrease at large |z| >> k, 
and for smaller H there is a substantial increase in magnitude 
and the lack of a flattening at large |z| >> h. Note that the last 
effects increase the contribution to g^ from the baryonic discs 
- perhaps dramatically - compared to a model that uses local 
disc properties only, and the second one mimics the effects of a 
constant DM density {g^{disc) const as |z| —> oo); both must 
result in a reduction in the amount of DM needed to fit the data. 
This is the opposite effect to that reported by Kuijken & Gilmore 
(1989b): in their Fig. 2, they assumed a smaller Rq = 7.4 kpc and 
a larger H = 4.5 kpc, i.e. H/Rq = 0.6. For values of H/Rq near 
unity, the gravitational field of the disc starts to approach that of 
a point mass, leading to the drop in g^(z) at large heights. 

The change in the form of g^ relative to that of an infinite 
disc is itself very similar to that of an infinite uniform disc with 
a different surface density and scale-height; in Fig. 1, the fits to 
the real form of g^(z) consists of that for an infinite disc with 
the nominal local surface density plus an additional infinite disc 
component with fitted surface density and scale-height (e.g. the 
decrease in g^(z) at large |z| for large N is then modelled as a 
negative additional surface density). When the Juric et al. (2008) 
thin- and thick-disc vertical scale-heights h of 300 and 900 pc 
and the corresponding radial scale-lengths N from Bensby et al. 
(2011) of 3.4 and 1.8 kpc, respectively, are used, the gravity due 
to the thin disc is nearly unaffected whereas the effective thick- 
disc surface density increases by 140%, resulting in an increase 
in the total effective surface density of 37%, and the additional 
component has a large effective vertical scale-length of 2.6 kpc. 
The changes in the form of g^(z) for the very thin gaseous disc of 
the MW is negligible because the radial scale-length is roughly 
twice that of the disc (Bigiel & Blitz 2012). 

Thus, the effects on the vertical gravity profile of non- 
uniform stellar structures in the Galactic neighbourhood of the 
Sun are quite significant and must be included in all verti¬ 
cal Jeans analyses. The naive separation of the gravity into an 
asymptotically constant disc and increasing DM component (e.g. 
the classic Kuijken & Gilmore “K” and “F” models) is not nec¬ 
essarily an adequate representation. Fortunately, even though the 
local vertical gravity is not that of an infinite disc with constant 
surface density, the effects of a realistic axisymmetric mass- 
distribution for the disc of the MW can be modelled using ad¬ 
ditional infinite disc components whose surface densities and 
scale-heights depend simply upon the true values. 


5. The vertical Jeans equation revisited 

The vertical Jeans equation can be easily derived from moments 
of the collisionless Boltzmann equation in cylindrical coordi¬ 
nates for a steady, axisymmetric stellar tracer component i with 
density pfR, z) (see BT08; Eq. 4.222b, p. 353). For a steady-state 
situation (F^. = 0 for k - R,z, i.e. no net flux of stars in R or z), 
one obtains the vertical, axisymmetric Jeans equation 

£ ^ ^ ^ 

This is an exact and full description of the connection between 
the local vertical gravity and these local stellar velocity disper¬ 
sions within the very modest assumptions made. Admittedly, 
there is evidence that the vertical structure of the MW’s thin disc 
is not stationary (e.g. Widrow et al. 2012, Carlin et al. 2014), but 
these are second-order effects that should not seriously affect the 
velocity dispersions, particularly at larger heights. Admittedly, 
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the gravitational field of the MW is not axisymmetric owing to 
factors like the central bar and the spiral arms, but the effects 
of the former are subtle (e.g. Bovy et al. 2014) and the non- 
axisymmetric effects of the latter are small due to local pitch an¬ 
gles of less than 20° (Vallee 2005; Levine, Blitz & Heiles 2006). 

Fortunately, the mean radial surface brightnesses of disc stars 
in most spiral galaxies are known to be exponential (Freeman 
1970) and the local vertical distributions of stars in the local 
MW are also observed to be exponential starting at quite modest 
heights above the plane of the disc (Juric et al. 2008). In partic¬ 
ular, this holds for mono-abundance tracers (Bovy, Rix & Hogg 
2012a; Bovy et al. 2012b; K13) so that one can confidently as¬ 
sume a doubly exponential density profile for any homogenous 
tracer population, and simply measure the radial scale-lengths 
Hi and vertical scale-heights /i, in order to characterise this most 
important property to an adequate accuracy : 

Pi(R, z) ^ PiiRo, . (3) 


Given the observed properties of p,, it is not unreasonable to 
make a similar assumption about the radial variation of the ve¬ 
locity dispersion: 

(4) 


A simple exponential disc with a constant disc thickness h re¬ 
quires cTjj oc S, i.e. H ~ H. More complicated dynamical mod¬ 
els of the MW suggest that is larger than H by ~20% (e.g. 
Sanders & Binney 2015), so assuming they are equal should not 
be a major source of error. 

To solve the vertical Jeans equation fully, one needs to know 
the radial and vertical behaviour of ctr^. This component gen¬ 
erally has a smaller effect than cr^ because the /i,- are generally 
much smaller than the Hi, so this contribution was totally left out 
in Z13. Several possible models for ctrj, have been proposed (see 
Amendt & Cuddeford 1991). For example, one can assume that 
the principle axes of the velocity ellipsoid are aligned in spheri¬ 
cal or cylindrical coordinates so that the correlation between the 
radial and vertical motions can be expressed via a tilt angle utut 
such that 


o-Rz ~ 2 tan(2a„7,) (ctr - 


(5) 


(see BT08, problem 4.32, p.391), where cr\ = crRR, cr\ = crzz- For 
the special case of Stackel potentials, it can be shown that a,,/, 
is aligned with the potential gradients and so should be inde¬ 
pendent of the tracer population used (e.g. Statler 1989 Binney 
2012). This form was used by Siebert et al. (2008) and Cassetti- 
Dinescu et al. (2011) to analyse abundance-undifferentiated data 
from the (southern) RAVE experiment close to the Galactic 
plane. Biidenbender et al. (2014) have measured amt for the 
abundance-differentiated dwarf G stars in the (northern) SEGUE 
sample and saw that amt increases in magnitude with increasing 
z: they fitted amt as a linear function in z but, given the errors, an 
equally good fit to the data is simply 


tM2am,{z)) = (-1.88 + 0.19)^ = c^. 

Kq Kq 

Both fits imply that the tilt angle is essentially aligned with the 
local spherical coordinate system. If an anisotropy parameter 
P{R,z) = is defined (e.g. Pasetto et al. 2012b), the as¬ 

sumptions p = const plus a tilt angle perfectly aligned with the 
Galactic center result in 

I (ctr - crfj (6) 


O'Rz 


1 -H ^ 
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Fig. 2. cr^ and cr^ data from Budenbender et al. (2014; B14) and Pasetto et al. (2012a,b); for the B14 SEGUE data, each colour represents a 
different [Fe/H] and [tr/Fe] mono-abundance tracer population (and was chosen to reproduce the corresponding plots in B14), whereas the Pasetto 
et al. data are for undifferentiated RAVE stars. 


(Kuijken & Gilmore 1989c), implying that the assumption c ^ 
const is not unreasonable for « j5R^. 

The SEGUE data from the B14 Table 1 are shown in Eig. 2, 
covering heights 0.5 < z < 2.5 kpc. Clearly, there is a tight 
connection between cr\ and cr^ that appears to be quite linear 
within a given tracer population and very similar for different 
tracer populationsIf one generically assumes 

tr®(/?,z)^fl,+h,tr®(/?,z) (7) 

(i.e. Pi = 1 jbi for a, = 0), one obtains a single equation express¬ 
ing ctr^ in terms of cr^, z, and two dispersion-model parameters 
T] and y, 

(t'rSR, z) ~ ^ Ci I (at 4- {bi - 1 )crW (R^ 2)) (8) 

= -(77,+7,o-®(f;,z))| 

where j] = —cajl and y = —c{b - l)/2 (the sign was chosen to fit 
the northern SEGUE data of B14). While this expression was de¬ 
rived using the tilt angle formalism, it constitutes a simple, phe¬ 
nomenological model for ctr^ in terms of the model parameters y 
and 7} and an assumed scaling in z that could cover other origins 
for the behaviour of this dispersion component. The relatively 
simple assumption that ctr^ is proportional to a linear function in 
cr^ whose constant term might be non-zero should not be con¬ 
strued as a potential physical problem: the fact that the function 
is non-zero for cr^ = 0 is irrelevant given that the dispersions 
are nearly isothermal for any given tracer population, and that 
what is needed is merely a functional connection between the 
two dispersions measures. 

The usefulness of this model for ctr. depends, of course, 
upon how well it fits the data, at least locally. In order to augment 

Actually, all of the data are consistent with the single assumption 
cr^ = (68.73 ± 1.1 km s“*)(Tj, but there is no obvious reason for this 
relation and it is not very useful for the analysis that follows. 


the B14 data with observations at lower heights, the thin-disc 
RAVE results from Pasetto et al. (2012a) were modelled using 
Eq.j^for cTjj and the epicyclic approximation for ctrr (Amendt 
& Cuddeford 1991). The fits to crkkiRo^z) were made using the 
python MCMC object emcee (Eoreman-Mackey et al. 2013) 
with the parameters Lk and the crk values at each height bin (0, 
+200, and ±400pc). A standard probability function InPoc 
and 30 MCMC walkers were used; after a burn-in of 2000 itera¬ 
tions each, the walkers were well converged and the probability 
density functions (PDE) were calculated from a total of 150000 
samples. While the values of cr^(RQ,z) and crRR{RQ,z) are rea¬ 
sonably well constrained, the final PDEs of the scale-lengths 
have long tails; if the radial scale-lengths are restricted to be 
with the reasonable range 0 < Lk < lOkpc, Lj. = l.Oj^^jkpc, 
and Lr = 2.8 ;*;“j kpc (68.2% bounds), consistent with being 
equal to a thin-disc density scale-length of ~ 3 kpc. When the 
RAVE dispersion points are added to Eig. 2, it appears that rj = 0, 
P const is a reasonable assumption for the thin disc, but the B14 
thick-disc data shows a much flatter relation. Thus, Eq.j^appears 
to be an adequate - if not excellent - local phenomenological 
model that is able to model the potentially different behaviours 
of different tracers. 

With these assumptions, one can then rewrite the local verti¬ 
cal Jeans equation as the ordinary differential equation 


d(T^^{Ra,z) _ 
dz 
where 


gz{Ro,z) + 


1 z 
hi |z| 


2T2 


o-«(Ro,z)- 


Ai = 


Ro 


RoHi 




2.2kpc(25kpc) (^82kpcj 


4r/ 


(9) 


( 10 ) 


r, 


- 1/2 


is a typical scale-length where the effects of cross-talk between 
ctr^ and cr^ appear, and Ail hi is a measure of the effects of the 
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(Trz term on for a given tracer; larger values imply less effect. 
The T]i term in Eq.j^is not proportional to cr^ and so behaves like 
a gravity term. Indeed, the proportionality to z makes it look like 
a constant apparent mass density p® defined by 


P 


(0 

n 


^nGRoHi 


l.SmMopc ^ 


2000 km^ s-2 


( 11 ) 

Ro Hi r‘ 

8.2kpc/ \3kpc/ ’ 


which represents a small but tracer-dependent and hence poten¬ 
tially noticeable correction to any fitted dark matter density for 
reasonable values of rp and //, . 

The equation for gziRo,z) used by Z13 to close the equa¬ 
tion assumes an infinite uniform disc plus a constant DM den¬ 
sity. In Section|^ we saw that a generalisation of their model 
- the sum of effective infinite discs with different effective sur¬ 
face densities and effective exponential vertical scale-heights 
Ij (the latter chosen to be easily distinguishable from the tracer 
scale-heights hj), plus a constant DM density Pdm{Ro, 0) - is a 
reasonable functional model; 


^ -IpDMiRo, 0)Z - ^ ^j{Ro)il - (12) 


This assumption also permits one to solve the vertical Jeans 
equation ( Eq.[9] l analytically when the simple model for cr^, 
(Eqs.[8]and|l l|l is used, 


InGh; 


M 


2 [pdm + pf) hi^iz-, hi. Ad + ^ 2,T(z; hi, T,-, Ij), (13) 

;=i 


where 


^(z;hi,Ai) 

r(z-,hi,Ai,lj) 

and where 


doff- 


/r? \h, 


4 

hi 







(14) 


-D 


I d ^ d kl \ 

\hi Ij 2Ail 


e-Wi 


D{x) = 




erfi(x) 


(15) 


is Dawson’s function (erfi(x) is the imaginary error function). 
This solution is only good for reasonable heights «2Aj/hi, but 
this is not a critical constraint on the applicability of the solution 
since d,//!, >> 1 for any reasonable value of y, and /r, (Eq. and 
the whole procedure breaks down (e.g. the assumptions of expo¬ 
nential density profiles and constant DM densities) if carried out 
too far from the disc plane. 

The behaviour of this solution relative to that of Z13 is gov¬ 
erned by the ratio Ai/hi, as is shown in Fig. 3. The effects of 
Ai are just apparent for Ail hi ~ 16, but are quite significant for 
smaller values. Thus, y, can significantly change the form of the 
solution. In the limit d, —> co, i.e. y, —> 0, this solution is the 
same as that used by Z13 and B14. However, for small Ai/hi, 
the simple asymptotic behaviour of the Z13 solution (flat disc 
and linear DM contributions to cr^) is modified so that both mass 
components produce increasing cr^ with height, making both the 
baryons and the DM more effective at explaining positive disper¬ 
sion gradients and thus principally reducing the total amount of 
DM needed. 
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Fig. 3. Top: the function hi. At, I/) representing the effect of a 

vertically exponential disc on cr^ for the ratio of vertical scale-heights 
lj/hi = 0.5 (red), 1 (blue), and 2 (black) for scaled ctr, parameters d,/!!, = 
4 (highest), 8, 16, and 32 (lowest). Bottom: the same for the function 
VTfx; hi, Ai) representing the effect of a constant density. 


6. Thin-disc versus thick-disc kinematics 

Biidenbender et al. (2014) showed that the crz(Ro,z) profiles of 
the thin- and thick-disc tracers in the SEGUE G dwarf dataset 
look very different (their Fig. 2); the former show a large gradi¬ 
ent of ~ 6kms“*kpc“^ at relatively low Galactic heights (|z| < 
1.3 kpc) and the latter half the gradient at much larger distances. 
As emphasised by B14, the simple model for cr^ used by Z13 is 
incapable of simultaneously explaining both the relative magni¬ 
tudes and the shapes of the two curves: one expects a rapid rise 
in cr^ at heights z < hdisc to a constant value due to the disc and 
a subsequent linear rise due to the constant DM density (Eq.[^. 
The difference between different tracers is a simple scaling with 
the tracer scale-heights /i, at large heights /r, > Ij, so the relative 
magnitudes are fixed by the relatively well-determined vertical 
density profiles; they cannot be changed at will, for example by 
modifying the relative amounts of baryonic and DM. Given that 
the full vertical Jeans equation can now be used and that one can 
model the local vertical gravity more realistically, is it now pos¬ 
sible to explain these differences, or are they simply due to tracer 
inhomogeneities? 

At first glance, the new model would appear to be able to 
solve the problem entirely: different values of y, (or, equiva¬ 
lently, of Ai) permit either a flat or an increasing asymptotic ctj. 
profile in the same disc but for different tracer populations, ex- 
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actly as seen by B14. The relative magnitudes are affected by 
different values of either y,- or rji, and by the difference due to 
relative values of the vertical scale-heights /t,. In contrast to the 
Z13 model, a flat profile does not necessarily require unrealis¬ 
tic amounts of baryons, since a small exponential radial scale- 
length can increase the effective surface density used to calculate 
the gravity by large amounts. If we find the right values for the 
kinematic and gravity model parameters, there would appear to 
be nothing in the way of determining a robust value for the DM 
density. 

As we saw in Section|^ a short radial scale-length H x 2 kpc 
for the thick-disc means that the importance and effective scale- 
height of this gravity component is much larger than previously 
considered, placing its effects in a region where an increase in 
cr^ was thought to be due to DM alone (Fig. 1). Since the grav¬ 
itational effects of the thick-disc should gradually increase with 
height, whereas those of the thin disc eventually saturate, there 
should be an additional gradient in cr^ for the thick-disc tracer 
data which samples the upper regions. This is, unfortunately, ex¬ 
actly the opposite of what is seen in B14’s two tracer samples, 
so the explanation is not simply a question of the gravity model. 


The main kinematic parameters are 77 , and y,-, expressible as 
the effective mass-density p® (Eq. 11 1 and the scale-length A, 
(Eq. 13 I. The value of A^ is proportional to Hijji, so the smaller 
the value, the more significant the effects. While B14 did not 
publish ctr^ or ctrr values for their two representative thin- and 
thick-disc samples, one can take the kinematic data in their Ta¬ 
ble 1 and split them into two kinematic groups, corresponding 
to the clearly thin discs ([a/Fe] < 0.2) and clearly thick discs 
([o'/Ee] > 0.28), i.e. the intermediate disc tracer with [Ee/H]= 
-0.35 and [a/Ee]=0.28 is excluded. This makes it possible to at¬ 
tempt to fit ctrz and hence attempt to explain the differences seen 
byB14. 


The first model (“A”) uses the tilt-angle formalism connect¬ 
ing ctr^ with both ctr and cr^ (Eq.[^. This entails fitting the 
model parameters a and b describing crR{cr^) (Eq.[7]l and then 
using these results to fit crR,{z, ctr, cr^) via the parameter c, from 
which y, 77 , and A can finally be derived for the B14 thin, thick, 
and combined datasets (model “Al”). The fits were again made 
using emcee (Foreman-Mackey et al. 2013) following a similar 
fitting process to that described in Section 2. The value for 

the thin dataset is large (3.1 vs. 1.1 for the thick dataset), but the 
parameter c representing the tilt-angle behaviour is clearly dif¬ 
ferent for the two populations (-1.4 + 0.2 vs. -2.5 + 0.3) unlike 
what one would expect for a gravitational potential representable 
by Stackel functions. The ratio Ajh is also different (13 + 3 vs. 
3.5 + 0.5), suggesting that the thick-disc cr^ should have a larger 
gradient than the thin disc, again exactly the opposite of what 
is seen in B14. Assuming that the tilt-angle behaviour is com¬ 
mon to both tracers, one can fit a common c parameter (model 
“A2”); this fit shows similar values for the model parameters and 
the same unexpected Ajh behaviour between the thin- and thick- 
disc tracers. Since there is no reason why the crR(crj.)-relation 
should be the same for all tracers, one can use an extreme model 
which avoids the use of 77 by setting a and 77 equal to zero; this 
model (“A3”) provides adequate explanations for crR{z, (Jz) when 
applied to the thin- and thick-disc data separately but results in 
very bad fits for the combined data (yj = 9.2). Despite this dras¬ 
tic modification of the kinematic model, the same unexpected 
Ajh behaviour is seen. 


The alternative phenomenological model (“Bl”) is to adopt 
the mathematical form of crRz<(Z,a-z), but to fit the two model 
parameters y and 77 directly, ignoring the question of the implied 


effect on <jr((Jz) within the tilt-angle assumption. The values of 
Trci/ model are lower simply because we are now fitting 

only the ctrz data rather than the higher quality ctr data as well, 
but the fitted values of 77 are not significant and the same Ajh 
behaviour is seen. Finally, one can fit a one-parameter model 
(“B2”) by assuming 77 s 0; these fits are just as good (yj is 1.2 or 
better) and ctrz is well fit, but one sees the same Ajh behaviour: 
y is 1.3 + 0.2 and 0.8 + 0.1, and Ajh is 8.9 + 0.8 and 3.2 + 0.2 for 
the thin- and thick-discs, respectively. 

Thus, one must conclude that the differences in the kinematic 
behaviours of the B14 thin- and thick-disc samples are not sim¬ 
ply due to the effects of the cross-dispersion term ctrz, at least 
as modelled using the tilt-angle paradigm. The fact that a real¬ 
istic gravity model should also result in the opposite behaviour 
seen by B14 reinforces this conclusion. If the effects seen in the 
SEGUE thin- and thick-disc data by B14 are not simply due to 
the inadequacies of the kinematic model, they must be due to 
the data. For example, the larger gradient in the B14 cr^ thin- 
disc data is mainly due to the outermost z point at 1.25 kpc: if 
this point is lowered, then the rise at lower heights can be due 
to the disc scale-heights alone. Alternatively, the effect might be 
due to the same behaviour seen in the MBCM data - a gradi¬ 
ent produced by mixing the different kinematic behaviours of 
different tracers. If one compares the B14 sampling with that 
of Bovy, Rix et al. (2012a; their Fig. 4), one sees that the verti¬ 
cal scale-heights of the thin-disc mono-abundance tracers in the 
B14 thin-disc samples vary from about 200 to 400 kpc: using 
the same simple calculation made for the BT12 analysis of the 
MBCM data in Section]^ one can produce an increase from 18 
to 25kms“' over the height range of 0.40 to 1.25 kpc (Fig. 3 in 
B14) with a local density ratio between the thinest and the least 
thin isothermal tracer component because the ratio of the relative 
contributions of the former over this range is roughly 10 % at the 
upper regions of the thin-disc data. For a reasonable mixture of 
densities, one then expects an artificial gradient even in the thin- 
disc cTj data due solely to population mixture. A similarly strong 
effect for the thick-disc data is not expected, since the thick-disc 
sample is dominated by the two most extreme thick-disc tracer 
populations which show roughly the same asymptotic values of 
cTz, but an artificial increase in the inhomogenous vertical disper¬ 
sion data must be present here as well. 

If there is no problem in the vertical gradients of cTj, mea¬ 
sured for truly homogenous data, is there one in the relative mag¬ 
nitudes, as suggested by B14 (the dashed solution in their left 
Fig. 2 diagram)? The B14 “no dark matter” solution used a very 
large baryonic surface density to explain the thick-disc data, but 
the model overpredicted the cr^ level for the thin disc by about 
20 % because the asymptotic ratio cr^(thick) /cr^(thin) is simply 
hrhickihthin for the Z13 model, independent of the relative contri¬ 
butions of baryonic discs and the DM halo. In the model for cr^ 
presented here, the ratio also depends upon the relative values of 
Ai/hi and 77 ,, i.e. in a complicated fashion that depends upon the 
full details of the models and the data used. A full discussion of 
these possibilities must be relegated to a subsequent paper, but 
it suffices to say that critical uncertainties in the vertical scale- 
heights and the obviously different values of y, are enough to 
reduce this discrepancy to a level where the model is principally 
able to explain the observations. Thus, the simplest explanation 
for the thin- and thick-disc kinematic behaviour seen by B14 is 
that the tracer data binned to represent only a thin and a thick 
disc is not homogeneous enough to permit the simple kinematic 
analysis used. 
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7. Conclusions 

I have shown that many previous attempts at estimating the local 
DM density have not been as robust as traditionally perceived. 
The usual baryonic corrections to the local dynamic density un¬ 
derestimate or ignore the uncertainties in the traditional gaseous 
densities and the traditional densities did not include the sub¬ 
stantial amounts of “dark” gas in the form of optically thick 
HI and CO-dark molecular gas now known to be present in the 
ISM. The non-local effects of Galactic structure increase the ef¬ 
fective dynamic density by another 25-35%, resulting in a total 
dynamic baryonic density of 120- 160mMopc~^, leaving prac¬ 
tically no room for measuring DM using mid-plane kinematics 
alone. The mixing of different tracer populations with different 
vertical scale-heights, radial scale-lengths, and dispersion mag¬ 
nitudes produces gradients in cr^(z) that can be misinterpreted as 
being due to DM whenever the effects of kinematic inhomogene¬ 
ity are not properly taken into account. The gravitational effects 
of doubly exponential discs are quite different in form and mag¬ 
nitude from the infinite discs often used, particularly for radial 
scale-lengths H <3 kpc. Fortunately, one can easily model this 
effect using additional infinite disc components with fitted sur¬ 
face densities and scale-heights. For a simple model of the MW 
thin and thick discs, the increase in the effective total gravitating 
surface density is quite large, ~40%. 

Assuming only axisymmetry, stationarity, and planar sym¬ 
metry, it is possible to derive a full analytic solution of the ver¬ 
tical Jeans equation for homogeneous tracers, a reasonably re¬ 
alistic local Milky Way disc plus a constant DM density that 
removes many of the approximations present in previous anal¬ 
yses. The dispersion cross-term o-r^ can be been modelled us¬ 
ing a phenomenologically justified expression motivated by, but 
not dependent upon, a tilt-angle assumption. The additional de¬ 
pendence on cTj in the vertical Jeans equation can induce tracer- 
dependent and detectable additional vertical dispersion gradients 
that can systematically reduce the amount of both baryonic and 
dark matter needed to explain the observations. 

The ctr^ model parameters were estimated using the SEQUE 
G dwarf data analysed by B12 and B14, either by assuming a 
tilt-angle model or by fitting a simple phenomenological model 
dependent on a linear function in cr^ and a factor of z/R- Despite 
the additional complexity of form available, it is not possible to 
explain the different behaviour of the “thin” and “thick” tracer 
data of B14 using kinematical effects alone; while the cross¬ 
dispersion effects should be large enough to be detected in the 
vertical profile of cr^(Ro,z), the changes are opposite to those 
needed to explain the different dispersion behaviours in the two 
B14 samples. The simplest explanation is that these two inho¬ 
mogeneous (non-mono-abundance) samples also show the tracer 
population mixture effects plaguing other analyses. 

Thus, an analysis of local stellar dispersion data to extract 
a robust local value of poM should in fact be possible, but it 
would require considerably more care and better data than has 
been used to date (e.g. GAIA). In a subsequent paper, I will ap¬ 
ply the techniques outlined herein to carefully chosen data in 
order to estimate the local DM density. 
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